In [1]:
#%matplotlib inline
# Import all the programs we want to use. If this gives an error, then you need to add these to your python path.
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pandas import DataFrame, Series # for convenience
import pims
import trackpy as tp
#import holopy as hp
import os
#import av
import scipy
import scipy.optimize as sco
import seaborn
from matplotlib.backends.backend_pdf import PdfPages
#%matplotlib notebook
from __future__ import division # this makes mathematical division work better
%pylab inline
In [2]:
scaling = 0.08431 #um/pixel
mpp = scaling
fps = 11.935
moviename = 'tracer+janus_3%_H2O2_5(red)2016-06-14'
os.chdir('C:\\Users\\Viva\\Desktop\\EPJ folder\\analysis')
tm = pd.read_pickle('filtered_data_with_drift_subtracted_tracer+janus_3%_H2O2_5(red)2016-06-14_pickled.pkl')
tracks = tm['particle'].astype(int).unique()
print size(tracks)
tm.head()
Out[2]:
In [3]:
plt.axis('equal')
ax = tp.plot_traj(tm, mpp = scaling, legend=False)
In [4]:
imsd = tp.imsd(tm, scaling, fps, max_lagtime=1000)
In [5]:
imsd.head()
Out[5]:
In [6]:
fig, ax = plt.subplots()
ax.plot(imsd.index, imsd, 'k-', alpha=0.1) # black lines, semitransparent
ax.set(ylabel=r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]',
xlabel='lag time $\Delta{}t$ [s]')
ylim(0,1)
fig.set_size_inches(5,3)
plt.title(moviename + '\nMSD')
Out[6]:
In [7]:
tm[tm['particle']==5233].head()
Out[7]:
In [8]:
print max(tm[tm['particle']==5233].x)
print min(tm[tm['particle']==5233].x)
In [9]:
## Filter out trajectories that have imsd 0.0.
list_of_long_trajectories = []
list_of_short_trajectories = []
for i in range(size(tracks)):
max_x = max(tm[tm['particle']==tracks[i]].x)
min_x = min(tm[tm['particle']==tracks[i]].x)
max_y = max(tm[tm['particle']==tracks[i]].y)
min_y = min(tm[tm['particle']==tracks[i]].y)
if isnan(max_x):
list_of_short_trajectories.append(tracks[i])
else:
if isnan(max_y):
list_of_short_trajectories.append(tracks[i])
else:
list_of_long_trajectories.append(tracks[i])
# remove each of the undesired trajectories
t_extended = tm.copy()
for i in list_of_short_trajectories:
t_extended = t_extended[t_extended.particle != i]
# now t_extended only consists of trajectories of particles that travel further than the threshold.
In [10]:
t_extended.head()
Out[10]:
In [11]:
imsd = tp.imsd(t_extended, scaling, fps, max_lagtime=1000)
In [12]:
fig, ax = plt.subplots()
ax.plot(imsd.index, imsd, 'k-', alpha=0.1) # black lines, semitransparent
ax.set(ylabel=r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]',
xlabel='lag time $\Delta{}t$ [s]')
ax.set_xscale('log')
ax.set_yscale('log')
fig.set_size_inches(3,3)
plt.title(moviename + '\nMSD')
Out[12]:
In [13]:
fig, ax = plt.subplots()
ax.plot(imsd.index, imsd, 'k-', alpha=0.1) # black lines, semitransparent
ax.set(ylabel=r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]',
xlabel='lag time $\Delta{}t$ [s]')
ax.set_xscale('log')
ax.set_yscale('log')
fig.set_size_inches(3,3)
plt.title(moviename + '\nMSD')
Out[13]:
In [14]:
emsd = tp.emsd(t_extended, scaling, fps,max_lagtime=60)
In [15]:
# I wonder whether these datapoints are normally distributed.
first = True
for i in np.logspace(0,log(len(imsd)-1)/log(10),num=6, dtype=int, base=10):
#print i
fig = figure(figsize=[5,.7])
# plot a histogram of the MSD values for a single lag time
imsd.iloc[i].hist(bins=30, grid = False, histtype='stepfilled', color='.35'
,label=r'$\Delta$t = {:.3} s'.format(imsd.index[i]))
try:
# plot a vertical line showing the ensemble mean square displacement for a single lag time
axvline(x=emsd.iloc[i], color='b', linewidth=1.5)
#,label='mean {:.2}'.format(emsd_zero_perox.iloc[i]))
except IndexError:
pass
if first:
plt.title(moviename + '\ndistribution of MSD values', loc='left')
first = False
plt.ylabel('Unweighted\noccurences')
plt.legend(frameon=True, fontsize='small', markerscale= 1)
plt.xlabel(r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]')
Out[15]:
In [16]:
res = tp.utils.fit_powerlaw(emsd) # performs linear best fit in log space, plots
res
Out[16]:
In [17]:
t0 = frange(.0666,10,.1)
fit = res.A[0]*(t0**res.n[0])
In [18]:
ax =emsd.plot(loglog=True, style='b.',linewidth=.3, grid=False, figsize = [3,3], label="emsd")
# plot fit
loglog(t0,fit, 'b')
ax.set(ylabel=r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]', xlabel='lag time $\Delta t$ [s]')
#ax.set_xticks([.1,1,10])
ax.get_xaxis().set_major_formatter(matplotlib.ticker.FormatStrFormatter('%g'))
plt.title(moviename + "\n2D MSD")
#savefig('/home/viva/group/viva/Analysis/2016-03-02/MSD_of_tracers_in_droplets_containing_both_tracers_and_Janus_particles.pdf')
#savefig('/home/viva/group/viva/Analysis/2016-03-02/MSD_of_tracers_in_droplets_containing_both_tracers_and_Janus_particles.png')
Out[18]:
In [19]:
# A new version of tp.motion.emsd() that calculates standard deviation.
# This function is copied from trackpy. (Please see the trackpy license.)
# I [Viva] added the calculation of biased weighted standard deviation.
def my_emsd(traj, mpp, fps, max_lagtime=100, detail=False, pos_columns=None):
"""Compute the ensemble mean squared displacements of many particles.
Parameters
----------
traj : DataFrame of trajectories of multiple particles, including
columns particle, frame, x, and y
mpp : microns per pixel
fps : frames per second
max_lagtime : intervals of frames out to which MSD is computed
Default: 100
detail : Set to True to include <x>, <y>, <x^2>, <y^2>. Returns
only <r^2> by default.
Returns
-------
Series[msd, index=t] or, if detail=True,
DataFrame([<x>, <y>, <x^2>, <y^2>, msd, N, lagt,
std_<x>, std_<y>, std_<x^2>, std_<y^2>,
std_msd],
index=frame)
Notes
-----
Input units are pixels and frames. Output units are microns and seconds.
"""
ids = []
msds = []
for pid, ptraj in traj.reset_index(drop=True).groupby('particle'):
msds.append(tp.motion.msd(ptraj, mpp, fps, max_lagtime, True, pos_columns))
ids.append(pid)
msds = pd.concat(msds, keys=ids, names=['particle', 'frame'])
results = msds.mul(msds['N'], axis=0).mean(level=1) # weighted average
results = results.div(msds['N'].mean(level=1), axis=0) # weights normalized
# Above, lagt is lumped in with the rest for simplicity and speed.
# Here, rebuild it from the frame index.
if not detail:
return results.set_index('lagt')['msd']
# Calculation of biased weighted standard deviation
numerator = ((msds.subtract(results))**2).mul(msds['N'], axis=0).sum(level=1)
denominator = msds['N'].sum(level=1) # without Bessel's correction
variance = numerator.div(denominator, axis=0)
variance = variance[['<x>', '<y>', '<x^2>','<y^2>','msd']]
std = np.sqrt(variance)
std.columns = 'std_' + std.columns
return results.join(std)
detailed_emsd = my_emsd(t_extended, scaling, fps, detail=True, max_lagtime=130)
In [20]:
plt.errorbar(detailed_emsd.lagt,
detailed_emsd.msd,
yerr = detailed_emsd.std_msd,
capthick=0,
alpha = 0.7,
linewidth=.2,
label="biased weighted standard deviation")
Out[20]:
values below 0 will not be plotted on loglog plot
In [21]:
plt.errorbar(detailed_emsd.lagt,
detailed_emsd.msd,
yerr = detailed_emsd.std_msd,
capthick=0,
alpha = 0.7,
linewidth=.2,
label="biased weighted standard deviation")
loglog(detailed_emsd.lagt, detailed_emsd.msd, 'b.', label="ensemble msd")
loglog(t0,fit, 'r', alpha=.7, label="power law fit") # plot fit
ax2 = plt.subplot(111)
ax2.set(ylabel=r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]', xlabel='lag time $\Delta t$ [s]')
ax2.get_xaxis().set_major_formatter(matplotlib.ticker.FormatStrFormatter('%g'))
plt.title(moviename + "\n2D MSD")
plt.legend(loc=2, fontsize='medium')
Out[21]:
In [22]:
fig, ax = plt.subplots()
ax.plot(imsd.index, imsd, 'k-', alpha=0.1, label="MSD of each particle")
#loglog(t0,fit, 'r', alpha=.7, label="power law fit") # plot fit
plt.errorbar(detailed_emsd.lagt, detailed_emsd.msd, yerr = detailed_emsd.std_msd, capthick=0,
linewidth=1.5, alpha=.3,
label="ensemble MSD with standard deviation")
ax.set(ylabel=r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]',
xlabel='lag time $\Delta{}t$ [s]')
ax.set_xscale('log')
ax.set_yscale('log')
ax.get_xaxis().set_major_formatter(matplotlib.ticker.FormatStrFormatter('%g'))
plt.title(moviename + '\n2D MSD, comparing imsd and emsd')
Out[22]:
In [23]:
# I wonder whether these datapoints are normally distributed.
first = True
for i in np.logspace(0,log(len(emsd)-1)/log(10),num=7, dtype=int, base=10):
#print i
fig = figure(figsize=[5,.7])
# plot a histogram of the MSD values for a single lag time
imsd.iloc[i].hist(bins=30, grid = False, histtype='stepfilled', color='.35',
label=r'$\Delta$t = {:.3} s'.format(imsd.index[i]))
# plot a blue vertical line showing the ensemble mean square displacement for a single lag time
axvline(x=detailed_emsd.iloc[i].msd, color='b', linewidth=1.5)
#,label='mean {:.2}'.format(emsd_zero_perox.iloc[i]))
# plot a red horizontal line showing the standard deviation
xmin = (detailed_emsd.iloc[i].msd)-(detailed_emsd.iloc[i].std_msd)
xmax = (detailed_emsd.iloc[i].msd)+(detailed_emsd.iloc[i].std_msd)
plot([xmin,xmax],[5,5], 'r')
if first:
plt.title(moviename + '\ndistribution of MSD values\n'+ '$\Delta$t = {:.3} s'.format(imsd.index[i]))
first = False
else:
plt.title('$\Delta$t = {:.3} s'.format(imsd.index[i]), fontdict={'fontsize':'medium'}, loc='center')
plt.ylabel('Unweighted\noccurences')
#plt.legend(frameon=True, fontsize='small')
plt.xlabel(r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]')
Out[23]:
In [24]:
from scipy.optimize import curve_fit
In [25]:
import scipy
scipy.__version__
# need at least version 14.0 of scipy.
Out[25]:
In [26]:
def powerlaw(t, A, n):
return (A * (t**n))
In [27]:
params, pcov = curve_fit(powerlaw, detailed_emsd.lagt, detailed_emsd.msd,
#p0=[res.A[0], res.n[0]],
sigma=detailed_emsd.std_msd,
absolute_sigma=True)
In [28]:
yfit = powerlaw(detailed_emsd.lagt,params[0],params[1])
#for y in yfit:
# if y < 0:
# yfit.replace(y, nan, inplace=True)
In [32]:
loglog(detailed_emsd.lagt,detailed_emsd.msd, 'k.')
loglog(detailed_emsd.lagt,yfit, 'm-')
plt.errorbar(detailed_emsd.lagt,
detailed_emsd.msd,
yerr = detailed_emsd.std_msd,
ecolor = 'k',
capthick=0,
alpha = 0.7,
linewidth=.3,
label="biased weighted standard deviation")
plt.xlabel('lag time $\Delta{}t$ [s]')
plt.ylabel(r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]')
#pylab.axes().set_aspect('equal', 'datalim')
plt.title(moviename + '\ncurve fit')
#print 'Fit to power law, MSD ~ A * $\Delta$t$^n$.'
#print 'Coefficient A = ' + str(params[0]) + ' ± ' + str(pcov[0,0])
#print 'Exponent n = ' + str(params[1]) + ' ± ' + str(pcov[1,1])
plt.text(0.1,10,
'Fit to power law, MSD = A * $\Delta$t$^n$.'
+ '\nExponent n = {:.3} $\pm$ {:.3f}.'.format(params[1],sqrt(pcov[1,1]))
+ '\nCoefficient A = {:.3} $\pm$ {:.3f}.'.format(params[0],sqrt(pcov[0,0])))
# The entries on the diagonal of the covariance matrix \Sigma are
# the variances of each element of the vector \mathbf{X}.
# source: https://en.wikipedia.org/wiki/Covariance_matrix#Generalization_of_the_variance
Out[32]:
In [34]:
plot(detailed_emsd.lagt,detailed_emsd.msd, 'k.')
plot(detailed_emsd.lagt,yfit, 'm-')
plt.errorbar(detailed_emsd.lagt,
detailed_emsd.msd,
yerr = detailed_emsd.std_msd,
ecolor = 'k',
capthick=0,
alpha = 0.7,
linewidth=.3,
label="biased weighted standard deviation")
plt.xlabel('lag time $\Delta{}t$ [s]')
plt.ylabel(r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]')
#pylab.axes().set_aspect('equal', 'datalim')
plt.title(moviename + '\ncurve fit')
#print 'Fit to power law, MSD ~ A * $\Delta$t$^n$.'
#print 'Coefficient A = ' + str(params[0]) + ' ± ' + str(pcov[0,0])
#print 'Exponent n = ' + str(params[1]) + ' ± ' + str(pcov[1,1])
plt.text(0.02,20,
'Fit to power law, MSD = A * $\Delta$t$^n$.'
+ '\nExponent n = {:.3} $\pm$ {:.3f}.'.format(params[1],sqrt(pcov[1,1]))
+ '\nCoefficient A = {:.3} $\pm$ {:.3f}.'.format(params[0],sqrt(pcov[0,0])))
# The entries on the diagonal of the covariance matrix \Sigma are
# the variances of each element of the vector \mathbf{X}.
# source: https://en.wikipedia.org/wiki/Covariance_matrix#Generalization_of_the_variance
Out[34]:
In [35]:
detailed_emsd.lagt.max()
Out[35]:
In [31]:
pcov
Out[31]:
In [32]:
imshow(abs(pcov), cmap="gray", interpolation="nearest", vmin=0)
plt.colorbar()
plt.title('Covariance matrix, absolute values')
Out[32]:
In [33]:
print moviename
print 'Coefficient A = ' + str(params[0]) + ' ± ' + str(np.sqrt(pcov[0,0]))
print 'Exponent n = ' + str(params[1]) + ' ± ' + str(np.sqrt(pcov[1,1]))
In [43]:
# Try linear fit
def linear(Dt,t):
return 4*Dt*t
Dt, Dt_pcov = curve_fit(linear, detailed_emsd.lagt, detailed_emsd.msd)
print Dt
print Dt_pcov
print "better fit next"
Dt, Dt_pcov = curve_fit(linear, detailed_emsd.lagt, detailed_emsd.msd,
sigma=detailed_emsd.std_msd,
absolute_sigma=True)
print Dt
print Dt_pcov
Dt = Dt[0]
Dt_pcov=Dt_pcov[0][0]
print 'Dt = ' + str(Dt) +' ± '+ str(Dt_pcov) + ' um^2/s'
linearfit=4*Dt*detailed_emsd.lagt
In [46]:
plot(detailed_emsd.lagt,detailed_emsd.msd, 'k.')
plot(detailed_emsd.lagt,linearfit, 'm-')
plt.errorbar(detailed_emsd.lagt,
detailed_emsd.msd,
yerr = detailed_emsd.std_msd,
ecolor = 'k',
capthick=0,
alpha = 0.7,
linewidth=.3,
label="biased weighted standard deviation")
plt.xlabel('lag time $\Delta{}t$ [s]')
plt.ylabel(r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]')
#pylab.axes().set_aspect('equal', 'datalim')
plt.title(moviename + '\ncurve fit')
#print 'Fit to power law, MSD ~ A * $\Delta$t$^n$.'
#print 'Coefficient A = ' + str(params[0]) + ' ± ' + str(pcov[0,0])
#print 'Exponent n = ' + str(params[1]) + ' ± ' + str(pcov[1,1])
plt.text(0.02,20,
'Fit to line MSD $= 4D_t \Delta t$ \n $D_t = $' +str(Dt) + ' $\pm$ ' + str(Dt_pcov) + ' $\mu m^2/s$')
# The entries on the diagonal of the covariance matrix \Sigma are
# the variances of each element of the vector \mathbf{X}.
# source: https://en.wikipedia.org/wiki/Covariance_matrix#Generalization_of_the_variance
Out[46]:
In [ ]: